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1.  Introduction 


ring  missile  flight  tests  at  White 


Sands  Missile  Range  (WSMR) ,  position  data  on  the  current  status  of 
the  missile  is  transmitted  20  times  per  second  from  radar  sites  to 
Univac  1108  computer.  Consecutive  pairs  of  such  data  are  averaged 
times  per  second  and  computations  for  plotting  displays  such  as  cur¬ 
rent  position,  range  verses  altitude,  or  impact  prediction  are  based 
upon  this  averaged  data. 


a 

10 


In  the  event  that  a  missile  veers  from  Its  planned  trajec¬ 
tory,  it  will  be  necessary  to  terminate  thrust  to  prevent  the  missile 
from  impacting  in  a  populated  area.  For  this  reason,  the  Range  Safety 
Officer  (RSO)  requires  that  for  each  computational  cycle  (10  per  sec¬ 
ond)  an  instantaneous  Impact  prediction  (IIP)  of  the  missile  be  com¬ 
puted.  This  point  is  the  intersection  of  the  missile  trajectory, 
should  thrust  be  terminated,  with  the  Clarke  Spheroid  (of  1866)  model 
of  the  Earth  at  an  altitude  of  4000  feet^^y^ 

If  the  effects  of  atmospheric  drag  are  neglected  then  a 
vacuum  IIP  can  be  rapidly  computed  using  Kepler's  central  force  equa¬ 
tions,  since  the  missile  trajectory  is  an  ellipse.  This  method  is 
applicable  for  a  variety  of  low  drag  vehicles  whose  exit  from  and  re¬ 
entry  to  the  Earth's  atmosphere  occur  at  high  angles.  Since  approxi¬ 
mately  one  millisecond  is  required  to  compute  a  vacuum  IIP,  it  is 
always  computed  at  a  10  per  second  rate.  Details  of  this  computation 
are  given  in  [1]. 


For  low  altitude  missiles,  whose  trajectory  is  signifi¬ 
cantly  affected  by  atmospheric  drag,  a  drag-corrected  IIP  calculation 
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is  also  required  by  the  RSO.  Since  the  ordinary  differential  equa¬ 
tions  describing  the  missile  trajectory  in  the  atmosphere  are  not 
solvable  in  closed  form,  as  they  were  in  the  case  of  an  elliptical 
trajectory  for  a  missile  in  a  vacuum,  these  equations  must  be  numer¬ 
ically  integrated  until  the  resulting  trajectory  pierces  the  Earth's 
surface.  The  numerical  method  currently  being  used  is  a  fourth-order 
Runge-Kutta  (RK)  method  with  approximately  ten  steps  per  trajectory 
length  or  a  second-order  RK  method  with  approximately  twenty  steps 
per  trajectory  length.  The  fourth-order  RK  method  is  used  whenever 
the  ballistic  coefficient  $  of  the  missile  is  greater  than  200  lbs/ 
ft2,  whereas  the  second-order  RK  method  is  used  whenever  g  is  less 
than  or  equal  to  200  lbs/ft2.  In  either  case,  approximately  forty 
evaluations  of  the  equations  of  motion  are  required  per  drag  IIP 
computation.  Since  the  Runge-Kutta  method  is  a  single  step  method, 
the  step  size  at  each  integration  step  is  independent  of  the  step 
size  at  all  previous  steps.  The  step  size  used  is  adjusted  at  each 
step  so  that  it  decreases  as  drag  increases,  in  order  to  obtain  a 
more  accurate  drag  IIP  in  the  same  number  of  integration  steps.  A 
derivation  of  the  equations  of  motion  and  the  method  of  solution  is 
described  in  [2],  whereas  the  step  size  adjustment  is  described  in 
[3]. 


2.  Statement  of  Problem.  Some  missions  at  WSMR  involve 
several  simultaneous  missiles.  A  drag  IIP  computation  is  required 
for  each  such  missile.  Experimental  runs  have  indicated  that  only 
one  drag  IIP  could  be  computed  at  a  ten  per  second  rate  and  at  most 
four  drag  IIPs  at  a  five  per  second  rate  [4].  Until  recently  drag 
I IPs  were  computed  at  a  five  per  second  rate,  with  a  linear  extra¬ 
polation  of  the  last  two  computed  drag  IIPs  being  used  to  approxi¬ 
mate  the  drag  IIP  at  the  next  intermediate  time.  Additional  mission 
requirements  have  necessitated  the  implementation  of  a  variable  rate 
for  computing  drag  IIPs,  ranging  from  ten  per  second  to  two  per  sec¬ 
ond  [4].  At  a  two  per  second  rate,  four  drag  IIP  approximations  via 
linear  extrapolation  are  required  for  each  computed  drag  IIP,  in 
order  to  output  IIPs  at  a  ten  per  second  rate.  Consequently,  the 
output  drag  IIPs  at  a  two  per  second  rate  are  not  as  smooth  and  ac¬ 
curate  as  they  are  at  a  five  or  ten  per  second  rate. 

The  purpose  of  this  paper  is  to  present  an  alternative 
method  of  computing  drag  IIPs  and  alternative  methods  of  obtaining 
approximations  to  the  drag  IIPs  at  intermediate  times,  in  order  to 
decrease  computation  time  and/or  increase  accuracy.  In  order  to  put 
this  paper  in  its  proper  perspective,  we  briefly  mention  some  pre¬ 
vious  investigations  toward  achieving  these  goals. 
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3.  Previous  Investigations.  In  1965  a  method  for  obtain¬ 
ing  a  drag-corrected  Kepler  IIP  for  Athena  missiles  was  investigated 

[5] .  Briefly  this  method  consisted  of  using  position  and  velocity 
data  from  a  previous  Athena  mission  to  calculate,  in  non-real  time, 
a  table  of  differences  between  Kepler  IIPs  and  drag-corrected  IIPs 
as  a  function  of  velocity.  During  the  next  Athena  mission  a  linear 
extrapolation  of  this  data  was  used  to  obtain  a  drag-corrected  Kepler 
IIP  in  real  time. 

In  1975  various  numerical  integration  methods  for  computing 
drag  IIPs  using  digital  and/or  analog  computers  were  investigated 

[6] ,  which  included  Adams-Moulton,  Milne-Hamming,  Euler,  a  variable 
order  Adams  method  called  D1FSUB  developed  by  C.  W.  Gear  [7,  8,  9], 
and  2nd,  3rd,  and  4th  order  Runge-Kutta  methods.  The  conclusions 
reached  were  that  the  RK  methods  were  better  than  the  other  methods, 
except  possibly  the  method  of  Gear.  Large  errors  were  associated 
with  the  analog  solutions. 

Subsequently,  alternative  methods  for  expressing  the  equa¬ 
tions  of  motion,  using  Encke's  method  [10,  pp.  29-35]  and  two  dif¬ 
ferent  versions  of  the  method  of  variation  of  parameters  [10,  pp. 
116-120  and  11] ,  were  investigated,  as  well  as  alternative  numerical 
integration  methods,  such  as  an  improved  variable  order,  variable 
step  Adams  method  [12],  a  rational  extrapolation  method  [13,  14,  15, 
16],  and  a  Gauss-Jackson  (E2)  method  [10],  The  conclusions  drawn 
were  that  the  method  of  variation  of  parameters  took  about  twice  as 
long  as  Cowell's  method  for  the  computation  of  the  same  drag  IIPs, 
whereas  Encke's  method  took  about  four  times  as  long  [17,  p.  131  and 
18,  p.  15].  Furthermore,  the  alternative  numerical  integration 
methods  investigated  offered  little  if  any  improvement  over  the  vari¬ 
able  step  RK  method  currently  being  used  [17,  p.  134], 

Recently,  an  "f  and  g  series"  impact  predictor  algorithm, 
which  is  based  upon  a  Taylor-series-in-time  representation  of  a  mis¬ 
sile's  position  and  velocity,  was  developed  for  which  a  "two  to  ten 
fold  reduction  in  computer  execution  time  for  satellite  orbits  and  a 
seven  to  ten  fold  reduction  in  execution  time  for  ICBM  trajectories" 
could  be  achieved  over  conventional  numerical  integration  [19,  p. 

59].  A  second  report  uses  the  f  and  g  series  technique  to  determine 
the  "geographical  distribution  of  debris  impact  coordinates  that 
would  result  if  the  missile  were  destroyed,"  [20,  p.  287].  This  tech¬ 
nique  is  an  extension  of  the  classical  f  and  g  series  used  in  the 
solution  of  Kepler's  equations  of  motion  [21,  pp.  107-111]. 
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Such  a  Taylor  series  type  of  solution  to  an  Initial  value 
problem  is  "generally  impractical  from  a  computational  point  of  view," 
[22,  p.  365].  In  fact  "...  the  necessity  of  calculating  the  higher 
derivatives  makes  Taylor's  algorithm  completely  unsuitable  on  high¬ 
speed  computers  for  general  integration  purposes,"  [23,  p.  330]. 
However,  in  "comparison  with  fourth-order  predictor-corrector  and 
Runge-Kutta  methods,  the  Taylor  series  method  can  achieve  an  appre¬ 
ciable  saving  in  computer  time,  often  by  a  factor  of  100,"  [24,  p. 
389].  In  view  of  these  statements,  the  equations  necessary  to  imple¬ 
ment  the  Taylor  series  method  are  derived  in  this  paper  in  order  to 
determine  whether  or  not  the  claims  made  for  satellite  orbits  and 
ICBM  trajectories  are  equally  valid  for  short  and  medium  range  tra¬ 
jectories,  such  as  those  experienced  by  the  missiles  tested  at  WSMR. 


4.  Equation  of  Motion.  The  vector  equation  of  motion  of 
a  missile  is 


(4.1) 


r  +  r  +  r ,  , 
u  g  d  * 


where  r  is  the  total  accelerationof  the  missile,  r^  is  the  unper¬ 
turbed  (Keplerian)  acceleration,  r  is  the  perturbative  acceleration 

© 


due  to  higher  order  harmonics  of  the  Earth's  gravitational  field,  and 
rd  is  the  perturbative  acceleration  due  to  the  Earth's  atmospheric 
drag. 


Missiles  which  are  launched  from  one  end  of  WSMR  and  which 
impact  at  the  other  end  do  not  travel  more  than  143  miles,  whereas 
missiles  which  are  launched  from  Green  River,  Utah  and  impact  on 
WSMR  do  not  travel  more  than  500  miles.  For  such  short  and  medium 

range  missiles  the  perturbative  acceleration  r  is  insignificant 
- 

compared  with  rd.  Therefore,  r^  -  C)  in  this  paper. 

The  expression  for  ry  is 

(4.2)  ry  -  (-u/r3)r  , 

_  _ IL 

where  r  ■  (x,  y,  z)  is  the  position  vector  of  the  missile,  r  -  (r*r)z 
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is  the  distance  of  the  missile  from  the  origin  of  an  inertial  rec¬ 
tangular  coordinate  system,  and  p  is  a  gravitational  constant.  If  r 
is  measured  in  feet  and  time  in  seconds,  then  p  *  1.406559714  X  106 
ft3/sec2. 

The  coordinate  system  used  is  right-handed,  with  origin  at 
the  Earth’s  center  (geocentric),  x-y  plane  in  the  equatorial  plane, 
z-axis  positive  through  the  South  Pole,  and  the  x-axis  aligned  at 
106°20'  west  longitude  at  the  instant  missile  position  and  velocity 
data  is  obtained.  Since  this  data  is  obtained  in  a  non-inertial 
(relative)  coordinate  system  which  rotates  with  the  Earth,  the  veloc- 
•  •  •  — 
ity  components  x  ,  y  ,  and  z  of  the  relative  velocity  vector  v 
*"  •  •  • 
must  be  converted  to  corresponding  components  x,  y,  and  z  of  the  in¬ 
ertial  velocity  vector  r  by 

(4.3)  x  «  xr  +  wy,  y  =  yr  -  tox,  z  -  zr  , 

where  to  ■  7.29211583  X  10~5  rad/sec  is  the  Earth's  angular  rate  of 

rotation. 

The  vector  acceleration  for  atmospheric  drag  is  given  by 

(4.4)  rd  *  -(p(h)vr/20)vr  , 

where  v^  is  the  velocity  of  the  missile  relative  to  the  Earth's  at¬ 
mosphere  in  ft/sec,  vf  -  (vr  •  vf)  ^2,  B  is  the  ballistic  coeffi¬ 
cient  of  the  missile  in  lbs/ft2,  and  p(h)  is  the  density  of  the  at¬ 
mosphere  in  slugs/ft3  at  the  current  missile  altitude  h.  The  Earth's 
atmosphere  is  assumed  to  rotate  with  the  Earth,  the  effects  due  to 
the  wind  are  neglected,  and  0  is  assumed  to  be  a  constant  for  each 
missile. 

If  r  denotes  the  velocity  of  the  missile  in  an  inertial 
coordinate  system  and  w  •  (0,  0,  -w) ,  then  the  relative  velocity  v^ 
is  given  by 

(4.  5)  vr  ■  r  -  w  X  r  . 
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The  atmospheric  density  is  approximated  by 

2 

(4.6)  p(h)  -  Ao  exp  (Ajh  +  A2h  )  , 

where  A0  *»  8.1283549  X  10“2,  Ax  -  -3.0319838  X  10~5,  and  A2  - 
-1.6214665  X  10-10  for  h  <  70,000  ft;  A0  -  1.7237156  X  l(f 1 ,  Ai  - 

-5.4682878  X  10_5,  and  A2  -  3.7544187  X  10"11  for  70,000  ft  <  h  < 
150,000  ft;  and  Aq  *  0  for  h  >  150,000  ft.  The  altitude  h  is  given 
by 

(4.7)  h  -  r  -  R  +  4000  , 

where  R  =  20929831.0  -  71303. 68411(z/r) 2  (cf.  [1,  p.  41])  is  the 
Earth  radius  of  the  Clarke  Spheroid  of  1866  at  the  same  latitude  as 
the  missile  but  at  the  4000  ft  altitude  of  WSMR. 

5.  Taylor  Series  Method.  The  Taylor  series  expressions 
the  missile's  position  r  =  r(t)  and  velocity  r  ■  r(t)  at  time  t  are 

(5.1)  7  *  Z±“0  (Ti/i!)70(1)  +  E:(t)  , 

(5.2)  7  -  (Ti/i!)r£i+1)  +  I2(t)  , 

where  t  «  t  -  tj  is  the  time  interval  or  step  size  for  re-initiali¬ 
zation  of  the  series  and  7^  -  7*^  (t0)  is  the  ith  time  derivative 
of  the  position  evaluated  at  epoch  tg.  The  error  vectors  Ei(t)  and 

E2(t),  due  to  the  truncation  of  these  Taylor  series  at  the  ro”^  term, 
are  given  by 


(5.3) 

Ei(t)  -  (t11*1/  (rH-1) !)  r(nfl)ai), 

rt 

o 

1  A 

Tpy 

H 

1  A 

(5.4) 

E2(t)  -  (Tn/n!)  7(n+1)(C2)  , 

<  _  < 

tQ  “  C2  “ 

Following  a  rule  of  thumb,  stated  by  Moore  ]25],  of  choosing  n  to  be 
approximately  equal  to  the  number  of  significant  decimal  digits  that 
can  be  carried  by  the  computer,  we  choose  n  to  have  a  maximum  value 
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of  6.  However,  for  small  values  of  the  ballistic  coefficient  B, 
smaller  values  of  n  with  smaller  step  sizes  t  generally  are  prefer¬ 
red. 


The  derivatives  r in  (5.1)  and  (5.2)  are  obtained  by 
**  • 

differentiating  ry  in  (4.2)  and  r^  in  (4.4)  with  respect  to  time  and 

v 

adding.  Hence  from  (4.1),  with  r^  ■  0, 

(5. 5)  dm  r/dtm  -  dm  r'  /dtm  +  dm  rjdtm. 

u  a 

In  the  development  of  the  f  and  g  series  algorithm  in  [19], 

the  expressions  for  the  second  and  higher  derivatives  of  r  failed 

•  ^ 

to  include  the  acceleration  r^  and  its  derivatives  [19,  p.  23,  (6)], 

resulting  in  the  omission  of  a  term  in  the  expression  for  d2r/dt2 
[20,  p.  301,  (15)],  and  the  omission  of  terms  in  all  higher  deriva¬ 
tives  of  r.  Therefore,  it  was  decided  to  determine  these  higher  deri- 

» 

vatives  of  ru  and  r^  independently  from  [19]  and  [20],  in  order  to 
correct  this  omission. 

Applying  Leibnitz'  Theorem  for  finding  the  nth  derivative 
of  the  product  of  two  functions  to  (4.2),  we  obtain 


in  which 

1d(l/r3)/dt  -  3r/r4  , 

d2(l/r3)/dt2  -  3r/r4  -  4rz/r5  , 

.. 

d3(l/r3)/dt3  «■  3r’/r4  -  36  rr/r5  +  60  r3/r6  , 
d4  (1/r3) /dt4  -  3riv/r4  -  48  irh5  -  36  r2/r5 
+  360  rzr/r6  -  360  r3/r7  . 
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Since 


(5.8) 


r  -  CE  •  T)1/2, 


we  easily  obtain 


r  =  (r  •  r)  7  r  . 

V  =  (-r  +  r*r  +  r*r)/r, 
r  =  (-3  r  r  +  3  r  •  r  +  r  •  r1  )  /  r  . 


The  corresponding  higher  derivatives  of  r ,  are  determined 

a 

by  applying  Leibnitz'  Theorem  to  (4.4),  to  obtain 


(5.9) 


where 


(5.10) 


in  26  AU / 


jit/  \  ,m-k 
d  (pvr)  d 


dk(pv  )  k  /.\  ,k  dk_iv 

~ r-  =  I  - £ 

i-0' 1 ' 


dtk  dt1*'1  ' 


From  (4.5),  we  have 


(5.11)  d^/dt*  -  d“  r/dt“  -  «  X  dm  r/dt®  . 

Since  vf  -  (vr  •  vr)  we  easily  obtain 


y  =  <?r  *  V  7  vr  ’ 


(5.12) 


(-v  +  v 


•  v  +  V  •  v  )  /  v  , 
r  r  r  r 


i»  _  jjL  _  tl* 

rr  =  (-3vr  vr  +  3vf  *  vr  +  vr  •  v  )  /  vr  , 

•  _  ..  2  ••  ••  * 


•  ••  ••  *-  21m  21  _  __  ___  „ 

(-4v  v  -  3v  +  3v  •  v  +4v  •  v  +v  •  v  )/v. 

rr  r  rr  rrrr  r 


The  higher  time  derivatives  of  p  in  (5.10)  are  obtained  by 
repeated  applications  of  the  chain  rule  of  differential  calculus  to 
p(h(t) )  in  (4.6).  We  first  note  from  (4.7)  that  the  time  derivatives 
of  the  Earth's  radius  R  are  small  compared  with  the  time  derivatives 

of  r  »  r.  Therefore,  we  have 
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(5.13)  d^h/dt111  =*  dmr/dtm  +  e  , 

m 

where  |ej  <<  ]dmr/dtm|. 
m 

Applying  the  chain  rule  to  p(h(t))  and  using  (5.13),  we  obtain 
dp/dt  =  r3p /3h  , 

|  d2p/dt2  =  r232p/3h2  +  r  3p/3h  , 

(5.14)  <  d3p/dt3  =  r333p/3h3  +  3  r  r  32p/3h2  +  r* 3p/3h  , 

/  d4p/dt4  =  r434pyjj^4  +  6  r^r  33p/3h3 

V  +  (3  r2  +  4  r  r*)32p/3h2  +  rlv  3p/3h  . 

Differentiating  (4.6),  we  have 
3p/3h  =  (Ai  +  2A2h)p  , 

32p/3h2  =  2A2p  +  (Ax  +  2A2h)3p/3h  , 

33p/ph3  -  4A23p/3h  +  (Aj  +  2A2h)  32p/3h2  , 

\^34p/3h4  =  6A232p/3h2  +  (Aj  +  2A2h)33p/3h3  . 

All  quantities  required  for  the  computation  of  the  deriva¬ 
tives  r^^,  i  =  1,  ....  6,  in  the  Taylor  series  (5.1)  and  (5.2)  can 
now  be  determined  from  the  preceding  equations.  To  implement  the 
Taylor  series  method,  a  value  of  n  is  choosen  between  2  and  6  with 
lower  values  corresponding  to  smaller  values  of  6.  The  step  size  t 
is  determined  such  that  the  maximum  number  of  steps  permitted  per 
trajectory  integration  is  not  exceeded,  and  such  that  the  local  trun¬ 
cation  errors  Ej  in  (5.3)  and  E2  in  (5. 4) < are  not  exceeded.  There¬ 
fore,  given  the  position  rg  and  velocity  rg  of  a  missile  at  epoch  tg, 

we  obtain  the  position  r  and  velocity  r  of  a  missile  at  time  t  from 
(5.I)  and  (5.2).  This  process  is  repeated  until  impact  time  T,  for 
which  r (T)  »  R(T) . 

Since  the  preceding  computations  were  made  in  an  Inertial 
coordinate  system,  the  true  impact  coordinates  Xj,  y^.,  and  Zj  can 


389 


*KUZANEK 


be  obtained  from  the  coordinates  x,j,,  y  ,  and  z^,  of  the  position  vec¬ 
tor  r(T)  at  impact  by  a  rotation  through  wT  radians  as  follows 

!Xj  =  Xj,  cos  uT  -  y sin  u)T  , 
yj  =  xT  sin  wT  +  yT  cos  wT  , 

2 1  =  *7* 

6.  Methods  for  Approximating  Drag  IIPs  between  Computa¬ 
tional  Cycles.  Ten  times  a  second,  or  once  every  100  milliseconds, 
all  data  on  a  missile,  such  as  its  position,  velocity,  and  drag  IIP, 
is  updated.  However,  the  computation  of  drag  IIPs  for  several  mis¬ 
siles  cannot  be  completed  in  fewer  than  100  milliseconds  even  using 
the  Taylor  series  method,  because  much  of  the  computer's  time  is 
spent  making  many  other  computations  during  each  100  millisecond 
time  period.  For  example,  if  drag  IIPs  can  be  computed  only  once 
every  500  milliseconds,  then  approximations  to  the  drag  IIPs  are 
required  four  times  per  500  millisecond  computational  cycle.  The 
current  method  for  obtaining  these  approximations  is  by  a  linear 
extrapolation  of  previously  computed  drag  IIPs,  considered  as  fun¬ 
ctions  of  range  time.  Unfortunately,  this  method  does  not  account 

for  new  values  of  the  missile's  initial  position  rg  and  velocity  rg 

at  these  intermediate  times  and  is  much  less  accurate  than  the  fol¬ 
lowing  improved  methods. 

The  first  improved  method  for  approximating  drag  IIPs  be¬ 
tween  computational  cycles  is  by  a  quadratic  extrapolation  of  the 
components  of  previously  computed  drag  IIPs,  considered  as  functions 
of  their  corresponding  components  of  the  vacuum  IIPs.  Since  vacuum 
IIPs  are  always  computed  every  100  milliseconds  anyway,  no  additional 
computer  time  is  required  for  their  use  in  this  method.  In  fact,  it 
requires  about  the  same  amount  of  computational  time  as  the  current 
linear  extrapolation  method  (one  millisecond),  yet  is  more  accurate 
by  an  average  factor  of  13.  Furthermore,  it  does  account  for  new 
intermediate  time  values  of  the  missile's  initial  position  rg  and 
• 

velocity  rg  since  the  vacuum  IIPs  at  these  times  are  functions  of 

rg  and  rg  .  However,  if  consecutive  pairs  of  components  of  vacuum 

IIPs  are  "close"  together,  then  this  quadratic  extrapolation  method 
may  yield  erroneous  drag  IIP  approximations.  Therefore,  the  follow¬ 
ing  method  which  avoids  this  problem  is  recommended. 


390 


*KUZANEK 


The  second  improved  method  for  approximating  drag  IIPs  be¬ 
tween  computational  cycles  is  by  a  quadratic  extrapolation  of  the 
differences  between  previously  computed  drag  IIPs  and  their  cor¬ 
responding  vacuum  IIPs,  considered  as  functions  of  range  time.  It 
too  requires  about  the  same  amount  of  computational  time  as  the  cur¬ 
rent  linear  extrapolation  method,  yet  is  more  accurate  by  an  average 
factor  of  6. 


In  general,  higher  order  extrapolations  may  yield  increas¬ 
ingly  worse  approximations  to  the  drag  IIPs  as  the  order  is  increased. 
This  would  be  especially  true  if  the  radar-determined  position  and 
velocity  of  the  missile  were  not  following  a  smooth  trajectory,  in 
which  case  an  accurate  drag  IIP  would  be  needed  most.  In  fact,  "If 
polynomial  extrapolation  must  be  done  with  poorly  behaved  functions, 
then  very  low  degree  extrapolation  is  usually  the  safest,  but  even 
this  should  be  carried  out  only  for  values  of  x  very  close  to  the 
tabulated  region,"  [26,  p.  58]. 

7.  Conclusions.  The  use  of  the  Taylor  series  method  re¬ 
sulted  in  equivalent  drag  IIPs  being  computed  in  two  to  ten  times 
less  time  than  by  the  currently  used  Runge-Kutta  method.  The  use  of 
the  method  of  quadratic  extrapolation  of  the  differences  between 
previously  computed  drag  IIPs  and  their  corresponding  vacuum  IIPs 
resulted  in  approximations  of  drag  IIPs  between  computational  cycles 
being  computed  over  six  times  more  accurately  than  by  the  currently 
used  linear  extrapolation  method,  with  about  the  same  amount  of  com¬ 
putational  time,  and  without  the  possibility  of  erroneous  approxi¬ 
mations. 
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